62. Multigrid and Multilevel Methods

Multigrid (MG) and Multilevel (ML) algorithms provide preconditioners with optimal condition numbers \(\kappa (C^{-1} A) = O(1)\), and optimal computational complexity \(O(N)\).

They can be seen as extension of the two level overlapping domain decomposition method to more levels.

• Ulrich Trottenberg, Cornelius W. Oosterlee, Anton Schuller: Multigrid, Academic Press, 2001

• Wolfgang Hackbusch: Multi-Grid Methods and Applications, Springer, 1985

In short, both methods are sub-space correction methods where the space splitting is defined by all basis functions together from a hierarchy of grids. Multilevel methods are additive Schwarz methods, while Multigrid methods are the multiplicative Schwarz methods. Both can be implemented efficiently by recursive algorithms.

We will present different theories for their analysis, one is based on sub-space decomposition using the ASM lemma, the other one uses the smoothing-and-approximation properties.

62.1. Multilevel preconditioner

We have a sequence of hierarchically refined meshes, which lead to a sequence of nested finite element spaces

\[ V_0 \subset V_1 \subset \ldots V_L \]

of dimension \(N_l = \operatorname{dim} V_l, l = 0 \ldots L\). We think of mesh-sizes \(h_l = 2^{-l}\), leading to finite element space dimensions \(2^{dl}\), where \(d\) is the spatial dimension.

We have prolongation matrices

\[ P_l \in {\mathbb R}^{N_l \times N_{l-1}}. \]

If \(v_{l-1}\) is a finite element function in \(V_{l-1}\) represented by the coefficient vector \(\underline v_{l-1}\). Then \(\underline v_l = P_l \underline v_{l-1}\) is the coefficient vector of the same function represented by the basis functions of \(V_l\).

If \(A_l\) and \(A_{l-1}\) are discretization matrices by a Galerkin method, then there holds

\[ A_{l-1} = P_l^T A_l P_l \]

Let \(D_l = \operatorname{diag} A_l\) be the Jacobi preconditioner (or some similar, cheap and local preconditioner).

2-level method

A 2-level preconditioner involving levels \(l-1\) and level \(l\) is

\[ C_{2L}^{-1} = D_l^{-1} + P_l A_{l-1}^{-1} P_l^T \]

By the ASM - Lemma we have the tools to analyze that this preconditioner has optimal condition number. However, a direct inversion of the matrix \(A_{l-1}\) is, up to a constant factor, as expensive as the inversion of \(A_l\). The multilevel method is to replace the inverses recursively by a multilevel preconditioner. On the coarsest level we invert the matrix \(A_0\):

\[\begin{eqnarray*} C_{ML,0}^{-1} & = & A_0^{-1} \\ C_{ML,l}^{-1} & = & D_l^{-1} + P_l C_{ML,l-1}^{-1} P_l^T \qquad \text{for} \; l = 1, \ldots L \end{eqnarray*}\]

from ngsolve import *
from ngsolve.la import EigenValues_Preconditioner

mesh = Mesh(unit_square.GenerateMesh(maxh=0.3))

fes = H1(mesh,order=1, dirichlet=".*", autoupdate=True)
u,v = fes.TnT()
a = BilinearForm(grad(u)*grad(v)*dx)

class MLPreconditioner(BaseMatrix):
    def __init__ (self, fes, level, mat, coarsepre):
        super().__init__()
        self.fes = fes
        self.level = level
        self.mat = mat
        self.coarsepre = coarsepre
        if level > 0:
            self.localpre = mat.CreateSmoother(fes.FreeDofs())
        else:
            self.localpre = mat.Inverse(fes.FreeDofs())
        
    def Mult (self, x, y):
        if self.level == 0:
            y.data = self.localpre * x
            return
        hx = x.CreateVector(copy=True)
        self.fes.Prolongation().Restrict(self.level, hx)
        cdofs = self.fes.Prolongation().LevelDofs(self.level-1)
        y[cdofs] = self.coarsepre * hx[cdofs] 
        self.fes.Prolongation().Prolongate(self.level, y)
        y += self.localpre * x

    def Shape (self):
        return self.localpre.shape
    def CreateVector (self, col):
        return self.localpre.CreateVector(col)

With operator notation we can define the multilevel preconditioner also as follows:

def MLPreconditioner2(fes, level, mat, coarsepre):
    prol = fes.Prolongation().Operator(level)
    localpre = mat.CreateSmoother(fes.FreeDofs())
    return localpre + prol @ coarsepre @ prol.T

a.Assemble()
pre = a.mat.Inverse(fes.FreeDofs())

for l in range(8):
    mesh.Refine()
    print ("ndof = ", fes.ndof)
    a.Assemble()
    pre = MLPreconditioner(fes,l+1, a.mat, pre)
    
    lam = EigenValues_Preconditioner(a.mat, pre)
    print ("lammin, lammax=", lam[0], lam[-1], \
            "kappa=", lam[-1]/lam[0])

ndof =  61
lammin, lammax= 0.48994111146454566 1.9701141796488941 kappa= 4.021124444445525
ndof =  217
lammin, lammax= 0.3922476241688647 3.3688433700649734 kappa= 8.588562842676819
ndof =  817
lammin, lammax= 0.38075965235975384 4.466009467399454 kappa= 11.72920880592628
ndof =  3169
lammin, lammax= 0.350072782680998 5.334434503847608 kappa= 15.238072674471765
ndof =  12481
lammin, lammax= 0.34830254682234674 6.024148612418195 kappa= 17.295735180170354
ndof =  49537
lammin, lammax= 0.34405057892289925 6.5755116439316215 kappa= 19.112049351921524

ndof =  197377

lammin, lammax= 0.3488640879128926 7.020027634327603 kappa= 20.1225287369803

ndof =  787969

lammin, lammax= 0.3524289346397199 7.381846752527455 kappa= 20.945631947256967

62.2. Multigrid Preconditioning

The multigrid preconditioner combines the local preconditioner and the coarse grid step sequentially. The multigrid preconditioning action is define as follows:

Multigrid \(C_{MG,l}^{-1} : d \mapsto w\)

• if l = 0, set \(w = A_0^{-1} d\) and return

• w = 0

• presmoothing, \(m_l\) steps:

  \[\hat w = w + D_{pre}^{-1} (d - A w)\]

• coasre grid correction:

  \[ \hat w = w + P_l C_{MG,l-1}^{-1} P_l^{T} (d - A w) \]

• postsmoothing, \(m_l\) steps:

  \[\hat w = w + D_{post}^{-1} (d - A w)\]

If the preconditioners \(D_{pre}\) and \(D_{post}\) from the pre-smoothing and post-smoothing iterations are transposed to each other, the overall preconditioner is symmetric. If the pre- and post-smoothing iterations are non-expansive, the overall preconditioner is positive definite. This can be obtained by applying forward Gauss-Seidel for pre-smoothing, and backward Gauss-Seidel for post-smoothing.

from ngsolve import *
from ngsolve.la import EigenValues_Preconditioner

mesh = Mesh(unit_square.GenerateMesh(maxh=0.3))

fes = H1(mesh,order=1, dirichlet=".*")
u,v = fes.TnT()
a = BilinearForm(grad(u)*grad(v)*dx)

class MGPreconditioner(BaseMatrix):
    def __init__ (self, fes, level, mat, coarsepre):
        super().__init__()
        self.fes = fes
        self.level = level
        self.mat = mat
        self.coarsepre = coarsepre
        if level > 0:
            self.localpre = mat.CreateSmoother(fes.FreeDofs())
        else:
            self.localpre = mat.Inverse(fes.FreeDofs())
        
    def Mult (self, d, w):
        # print ("level = ", self.level)
        if self.level == 0:
            w.data = self.localpre * d
            return
        
        prol = self.fes.Prolongation().Operator(self.level)

        w[:] = 0
        self.localpre.Smooth(w,d)
        res  = d - self.mat * w
        w += prol @ self.coarsepre @ prol.T * res
        self.localpre.SmoothBack(w,d)


    def Shape (self):
        return self.localpre.shape
    def CreateVector (self, col):
        return self.localpre.CreateVector(col)

a.Assemble()
pre = MGPreconditioner(fes, 0, a.mat, None)

for l in range(5):
    mesh.Refine()
    fes.Update()
    print ("ndof = ", fes.ndof)
    a.Assemble()
    pre = MGPreconditioner(fes,l+1, a.mat, pre)
    
    lam = EigenValues_Preconditioner(a.mat, pre)
    print ("lammin, lammax=", lam[0], lam[-1], \
            "kappa=", lam[-1]/lam[0])

ndof =  61
lammin, lammax= 0.7546521555146197 0.9984968279372552 kappa= 1.3231219451779748
ndof =  217
lammin, lammax= 0.5814333684266482 0.996629519645289 kappa= 1.7140906830685634
ndof =  817
lammin, lammax= 0.5486872840893804 0.9969908648358896 kappa= 1.817047512756795
ndof =  3169
lammin, lammax= 0.4985451720975005 0.9969259683462246 kappa= 1.9996702889567963
ndof =  12481
lammin, lammax= 0.4886987998302721 0.9964727026292579 kappa= 2.0390324326053975

f = LinearForm(1*v*dx).Assemble()
gfu = GridFunction(fes)
from ngsolve.krylovspace import CGSolver
inv = CGSolver(mat=a.mat, pre=pre, printrates=True)
gfu.vec.data = inv * f.vec

CG iteration 1, residual = 0.1748470807820223     
CG iteration 2, residual = 0.010391750485190051     
CG iteration 3, residual = 0.0013935264442234666     
CG iteration 4, residual = 0.0002485494511032407     
CG iteration 5, residual = 5.072626843296547e-05     
CG iteration 6, residual = 9.225744503561945e-06     
CG iteration 7, residual = 1.4178791475215712e-06     
CG iteration 8, residual = 2.3661511457730928e-07     
CG iteration 9, residual = 4.504902220776862e-08     
CG iteration 10, residual = 7.857866119833598e-09     
CG iteration 11, residual = 1.5078913529860317e-09     
CG iteration 12, residual = 2.937151908609418e-10     
CG iteration 13, residual = 4.8061282523741893e-11     
CG iteration 14, residual = 9.233582022146966e-12     
CG iteration 15, residual = 1.8488623475459732e-12     
CG iteration 16, residual = 3.510233606690287e-13     
CG iteration 17, residual = 7.073751997487315e-14     

from ngsolve.webgui import Draw
if fes.ndof < 20000:
    Draw(gfu, order=1)

62.3. Projection matrices from the finest level

It is often not feasible to assemble matrices on the coarse level, for example when solving non-linear problems. Then it is useful to calculate coarse grid matrices from the matrix on the finest level using the Galerkin property

\[ A_{l-1} = P_{l}^T A_l P_l \]

[ Requires NGSolve updates from 14.04.2021, still tricky: freedofs on coarser levels]

from netgen.geom2d import unit_square
from ngsolve import *
from ngsolve.la import EigenValues_Preconditioner

mesh = Mesh(unit_square.GenerateMesh(maxh=0.3))

fes = H1(mesh,order=1, dirichlet=".*", autoupdate=True)
u,v = fes.TnT()
a = BilinearForm(grad(u)*grad(v)*dx)

for l in range(5):
     mesh.Refine()
    
a.Assemble();

class ProjectedMG(BaseMatrix):
    def __init__ (self, fes, mat, level):
        super(ProjectedMG, self).__init__()
        self.fes = fes
        self.level = level
        self.mat = mat
        if level > 0:
            self.prol = fes.Prolongation().CreateMatrix(level)
            self.rest = self.prol.CreateTranspose()
            coarsemat = self.rest @ mat @ self.prol # multiply matrices
            self.localpre = mat.CreateSmoother(fes.FreeDofs())
                
            self.coarsepre = ProjectedMG(fes, coarsemat, level-1)
        else:
            self.localpre = mat.Inverse(fes.FreeDofs())
        
    def Mult (self, d, w):
        if self.level == 0:
            w.data = self.localpre * d
            return
        w[:] = 0
        self.localpre.Smooth(w,d)
        res = d - self.mat * w
        w += self.prol @ self.coarsepre @ self.rest * res
        self.localpre.SmoothBack(w,d)

        
    def Shape (self):
        return self.localpre.shape
    def CreateVector (self, col):
        return self.localpre.CreateVector(col)

# print (fes.mesh.levels)
pre = ProjectedMG(fes, a.mat, fes.mesh.levels-1)

print ("ndof: ", fes.ndof)
lam = EigenValues_Preconditioner(a.mat, pre)
print ("lammin, lammax=", lam[0], lam[-1], \
            "kappa=", lam[-1]/lam[0])

ndof:  12481
lammin, lammax= 0.48769694885003023 0.9964653220989623 kappa= 2.043205979550595

f = LinearForm(1*v*dx).Assemble()
gfu = GridFunction(fes)
from ngsolve.krylovspace import CGSolver
inv = CGSolver(mat=a.mat, pre=pre, printrates=True)
gfu.vec.data = inv * f.vec

from ngsolve.webgui import Draw
if fes.ndof < 20000:
    Draw(gfu, order=1)

CG iteration 1, residual = 0.1748470807820223     
CG iteration 2, residual = 0.010391750485190443     
CG iteration 3, residual = 0.0013935264442234423     
CG iteration 4, residual = 0.0002485494511032409     
CG iteration 5, residual = 5.0726268432965385e-05     
CG iteration 6, residual = 9.225744503561916e-06     
CG iteration 7, residual = 1.4178791475215523e-06     
CG iteration 8, residual = 2.366151145773084e-07     
CG iteration 9, residual = 4.5049022207766795e-08     
CG iteration 10, residual = 7.857866119833164e-09     
CG iteration 11, residual = 1.5078913529859374e-09     
CG iteration 12, residual = 2.937151908609169e-10     
CG iteration 13, residual = 4.806128252373669e-11     
CG iteration 14, residual = 9.233582022145108e-12     
CG iteration 15, residual = 1.848862347545574e-12     
CG iteration 16, residual = 3.510233606689523e-13     
CG iteration 17, residual = 7.073751997485764e-14